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Abstract. We present similarity solutions for adiabatic 
bubbles that are blown by winds having time dependent 
mechanical luminosities and that are each mass-loaded at 
a rate per unit volume proportional to T 5 / 2 r x , where T is 
the temperature, r is the distance from the bubble center, 
and A is a constant. In the limit of little mass loading a 
similarity solution found by Dyson ( 1973 ) for expansion 
into a smooth ambient medium is recovered. For solutions 
that give the mass of swept-up ambient gas to be less than 
the sums of the masses of the wind and the evaporated ma- 
terial, A > 4. The Mach number in a shocked mass-loaded 
wind shows a radial dependence that varies qualitatively 
from solution to solution. In some cases it is everywhere 
less than unity in the frame of the clumps being evapo- 
rated, while in others it is everywhere greater than unity. 
In some solutions the mass-loaded shocked wind under- 
goes one or two sonic transitions in the clump frame. Max- 
imum possible values of the ratio of evaporated mass to 
stellar wind mass are found as a consequence of the evapo- 
ration rates dependence on temperature and the lowering 
of the temperature by mass-loading. Mass-loading tends to 
reduce the emissivity in the interior of the bubble relative 
to its limb, whilst simultaneously increasing the central 
temperature relative to the limb temperature. 



Observations of wind-blown bubbles (WBBs) (e.g. 



Smith et al. 1984; Meaburn et al. 1991) have led to the 
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conclusion that the interactions of the winds with clumps 
of stellar material ejected during prior stages of mass- 
loss greatly affect the structure and evolu tion of the bub- 
bles (see Williams, Hartquist & Dyson 1995 and refer- 
ences therein). Consequently, a number of numerical in- 
vestigations of this process have been conducted. De- 
tailed one-dimensional, time dependent hydrodynamic 
models of specific W BBs have bee n constructed by Arthur, 
Dyson fc Har tquist ( |1993| , [I994D and Arthur, Henney & 
Dyson (1996). Numerical solutions for the Mach num- 
bers of steady mass-loaded winds have been presented 
by Williams et al. ( 1995| ), and studies of steady hydro- 
magnetic mass-loaded winds have been performed by 
Williams, Dyson, & Hartquist (|1999D . 

In this paper we derive similarity solutions for the 
structures and evolution of mass-loaded WBBs. This work 
complements that on similarity solutions for SNRs ob- 
tained by Chieze & Lazareff ( |1981| ) and Dyson & Hartquist 
(1987). Here we consider a constant rate of energy input, 
E, rather than a constant total energy, and assume that 
the mass-loading occurs due to conductively driven evap- 
oration. The solutions are potentially relevant to WBBs 
created by a fast wind interacting with a clumpy AGB 
superwind, by the wind of a young high-mass star inter- 
acting with surrounding molecular material, and the wind 
of an active galactic nucleus impacting its environment. 



1. Introduction 

Many authors have pointed out that mass pick-up from 
embedded clumps can significantly modify the structures 
of a wide variety of astrophys ical fl ows. These include, for 
example, McKee & Ostrikcr ( 1977 ) i n the ir theory of the 
interstellar medium, Hartquist et al. ( 1986 ) i n the study of 
wind blown nebulae, and Chieze & Lazareff (1981), Cowie, 
McKee & Ostriker ( |1981| ) and Dyson & Hartquist ( [198 7ft 
in the study of supernova remnants (SNRs). 



2. Similarity Solutions 

2.1. Overview of assumptions 

The evaporation of spherical clouds by hot gas has been 



analyzed by Cowie & McKee (1977) and McKee & Cowie 
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(1977). If the mean free path for electron-electron energy 
change in the ambient gas is less than about half of the 
radius of a clump, the mass evaporation rate from an in- 
dividual clump scales as T 5 / 2 , where T is the tempera- 
ture of the material surrounding the clump. For a typical 
jump of 1000 km s -1 at the reverse shock (corresponding 
to a post-shock temperature of T w 1.4 x 10 7 K), and 
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a clump electron number density of n e w 10 4 , the mean 
free path for electron energy exchange is A ss 5 x 10 15 cm. 
This is somewhat smaller than the lengthscale of a clump 
(e.g. in a planetary nebula, Meaburn et al. 1998), so the 



assumption of unsaturated conduction with a T 5 / 2 tem- 
perature dependence of p, the mass evaporation rate per 
unit volume, should be reasonable throughout most of the 
mass- loading region if mass- loading is significant. 

We also assume that the clumps are stationary with 
respect to the star, and that they are also cold and dense. 
Therefore the clumps have neither kinetic or internal en- 
ergy, and the picked-up mass is thus injected at zero ve- 
locity. The stripped mass then acquires momentum and 
energy from the general more tenuous flow. Hence the only 
source term in the hydrodynamic equations is p. In addi- 
tion we assume that at any given point the injected mass 
reaches the general flow velocity and temperature instan- 
taneously i.e. that the characteristic scale length of the 
mixing region is much smaller than the dimensions of the 
bubble. We also assume that the clumps are sufficiently 
numerous that they can be considered continuously dis- 
tributed such that a hydrodynamic treatment applies. We 
suppose that conductive energy transport can be neglected 
in the energy equation describing the global evolution of 
the WBB. The WBBs are assumed to evolve adiabati- 
cally, whilst the swept-up interclump medium is assumed 
to catastrophically cool into a thin dense shell. The clumps 
are assumed to pass through this shell without any in- 
teraction, and then on into the region of shocked stellar 
wind. The unshocked stellar wind is assumed to have a 
flow speed much greater than its thermal speed. Hence its 
thermal energy is negligible with respect to its kinetic en- 
ergy, and since no mass-loading occurs in this region, we 
ignore the energy equation for this part of the flow. Fi- 
nally we assume that any external forces, such as gravity 
or radiation pressure, are negligible. 

A schematic of the flow topology is shown in Fig. |l|. 
The hydrodynamical equations for the SSW region are: 
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where the symbols have their usual meanings. For con- 
ductively induced mass-loading, the evaporation of the 
clumps only occurs in the SSW region. In the USW re- 
gion we assume no interaction between the clumps and 
the fast stellar wind. 

2.2. The Similarity Variables 

Let the interclump medium have a density of the form 
pic = Po T ® : and let us consider the case where the mass- 
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Fig. 1. Schematic of the flow structure in our solutions. 
The central star blows out a bubble into the surround- 
ing medium which is characterized by a region of un- 
shocked stellar wind (USW), a region of shocked stel- 
lar wind (SSW), and a dense shell of swept-up material 
(SUM). We assume that the shocked ambient medium 
quickly radiates all of its internal energy, whilst the cool- 
ing time of the SSW is assumed to be much longer than 
the current age of the bubble. Mass-loading only occurs 
in the SSW. 

loading rate is also radially dependent: p = Q (e/p) 5 ^ 2 r x . 
The problem as its stands has five parameters with inde- 
pendent dimensions, namely r, t, Q, po, and the rate of 
energy input by the stellar wind E. The solution inside 
the bubble depends on r, t, Q, and E, whilst the solu- 
tion ahead of the contact discontinuity depends on r, t, 
po, and E. Inside the bubble we can construct a variable 
r = x E u t v Q w with a dimension of length. Outside the 
bubble we construct a similar variable r' = x' E x t v pg. 
Both x and x' are dimensionless. For a similarity solution 
to exist v = y, which demands that A = (5 + 70) /3. This 
reduces the number of dimensional parameters to 4: r, t, 
Q, and E. The solution may then be expressed in terms 
of the dimensionless similarity variables x, f(x), g(x), and 
h(x) where: 

r = x Q-VC10+A) £?V(10+A) t 7/(10+A) (4) 
p = Q 5 /(10+A) £>(5+A)/(10+A) t -(5-3A)/(10+A) (g) 
U = Q-V(10+A) £l/(10+A) t -(3+A)/(10+A) g ^ (g) 

e = QV^+x) £(7+a)/ ( io+a) 4 -(ii_a)/(io+a) h{x) (7) 
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The hydro-dynamic equations for the region of shocked 
stellar wind become the following set of coupled ordinary 
differential equations: 



fi and gi , the Rankinc-Hugoniot conditions for a 7 = 4/3 
gas specify that the post-shock value of the thermal energy 
density, h,2, is: 



x 10 + A 

h\ 5/2 

f) 

3 + A 



x x = 



-9 



10 + A 
g /fcN 5/2 
/ \f 



x x = 



(8) 



(9) 



(-9 - Wl) 



10 + A 



10, g 2 fh^ 5/2 



3x 



2 V/ 



x x = (10) 



where a prime denotes derivation with respect to x. It is 
simple to rearrange these equations to find /', g', and b! . 

2.2.1. Boundary conditions 

We start our integration at a point very close to the star. 
As r — > 0, 4nr 2 pv — > M c , where M c is the current rate of 
mass- loss of the star {M is a function of time if A ^ —3). 
Substituting the similarity variables we obtain: 



M c = lim 4irx 2 fgZ 

x—>0 



(11) 



where E = Q 2 /(io+a) £>(8+a)/(io+a) t (6+2\)/(w+\) _ We de _ 
fine </>, such that 



M, 



lim 47r a; /<? 



(12) 



which is a measure of the ratio of evaporated mass to wind 
mass. We also require that as r — > 0, 2irr 2 pv 3 — > E. i.e. 
that 



£ = lim 2ttx 2 fg 3 E 

x—>0 



(13) 



Cancelling the E's and substituting into Eq. [12] we obtain 
the values of / and g at x = Ax, the inner boundary of 
our integration: g = y / 2/(f>, f = l/(27r Ax 2 <? 3 ). 

The inner shock is positioned at x — x- ls . Once the inte- 
gration has proceeded to this point the standard Rankine- 
Hugoniot jump conditions across a discontinuity £ must 
be satisfied: 

Hs=0; [pt' 2 +p] E = 0; [1^ + ^-^ = 0(14) 

2 J-l p 

The pre-shock flow has an effectively infinite Mach number 
and the strong shock jump conditions apply. If the post- 
shock values for the density and the velocity are given by 



s) 2 



(15) 



where g s = 7x- ls / (10 + A) is the velocity of the inner shock 
relative to the star. 

Once the post-shock conditions have been calculated, 
integration proceeds until the contact discontinuity (CD) 
is reached. This occurs once the flow velocity is equal to 
the coordinate velocity (i.e. v = dR/dt), at which point 
g(x = x c d) = 7cc/(10 + A) is satisfied. At the CD, we must 
also satisfy conservation of momentum. This requires that 
the pressure inside the bubble is equal to the impulse on 
the swept up shell: 



. 2 d(M sh v) dM sh ,-,<!« 
47rr p = __ =v __+ Msh - 



(16) 



where here M s h is the mass of the swept-up shell (M s h 
34 



47r Po^d P /(3+/3)), and dM sh /dt = 4irr 2 pv. Eq. || reduces 



to 



a 2 , Po 
por p v H 



r i+B dv 



3 + P dt 



= P 



(17) 



We now construct a dimensionless constant, 9, from Q, E, 
and po: 



Q-1/(10+A)£il/(10+A) 
p -l/{5+8)^ 1/{5+0) 



(18) 



It provides a measure of the ratio of the evaporated mass 
to the swept-up mass. Inserting the similarity variables 
into Eq. [l?] and substituting Eq. |l8|, we find after some 
rearrangement : 



(7 - l)h 



l/(5+/3) 



3+A 



(19) 



(10+A)(3+/3) 



xg 



The value of 6 is evaluated once the integration stops at 
the CD. 



2.3. Scale Transformation and Normalization 



The dimensionless Eqs. |8J-[L0| (and their counterparts for 
the USW region) are invariant under the following trans- 
formation, which we shall call a normalization (cf. Chieze 



Lazareff 1981) 
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These can be combined to obtain the normalizations for 
the mass, and for the kinetic and thermal energies, of the 
bubble: 



m = 47r / fx 2 dx 
k = 47r j — fg 2 x 2 dx 
= 4-7T / hx 2 dx 



,,8+Ai 



,10+Ai 



, 10+Ai 



(21) 
(22) 
(23) 



The power of a in the square brackets indicates how each 
integral scales in a normalization. The full equations for 
the mass, kinetic and thermal energies, in the bubble are 



4na 8+x Zt 



fx 2 dx = a s+x ^-m 



KE h = 2ira 



x Et / fg 2 x 2 dx = a w+X 

x=0 



Etk 



„!0+A T?. 



10 + A 77 



Eti 



(24) 
(25) 
(26) 



a:=0 



The mass in the USW (SSW) region, M usw (M ssw ), is 
obtained through the replacement of m in Eq. |i| with 
Tn^sv, (m ssw ), where the integral is evaluated from x = 
to x — Xi s (from x — X[ s to x = x cd ). The total mass 
carried from the star by its wind is 



M(t) dt 



10 + A . 

-Mr t 



16 + 3A 



(27) 



The degree of mass- loading in the bubble, <&b, can be mea- 
sured by the ratio Mb/M w . Starting with the expression 
for M s h, (cf. text following Eq. [If]), using Eq. |^ to express 
r cd in terms of x cd , employing Eq. [l^ to rewrite po in 
terms of 8, and using Eq. |Il] and the definition of S, we 
find 



M sh = 



— a s + x 6^x s +^ 



3 + /3 



(28) 



The ratio of the swept-up mass in the shell to the bubble 
mass is hence: 



M sh M h = - cd 
(3 + p) m 



(29) 



The kinetic energy of the shell is KE^ = iM s h^ d . Using 
Eqs. |2|, pi, M, and | through |o| we obtain 



KE * = ^-^ x Ete^xir 9 2 cd 



(30) 



We can also calculate the pdV work done on the shocked 
interclump medium in compressing it into a negligibly thin 
shell. The initial energy per unit mass of the interclump 



gas immediately after passing through the forward shock 
is: 



Eh. 



1 



3pi 

2/91 



(31) 



where here vi is the post-shock flow velocity (= jv cd ) and 
Pi and pi are the post-shock pressure and density. If we 
assume that the pressure and velocity of the surround- 
ing medium are both negligible, the post-shock pressure 
is pi = fpo^fd^cd- Substituting for p\ into Eq. [n], we find 
that the sum of the kinetic and thermal energies per unit 
mass behind a strong shock is E\ — j§v 2 A . From our ear- 
lier assumptions, the swept up gas ends up with a veloc- 
ity v — v cd (i.e. it is accelerated) and no internal energy. 
Hence, the final energy per unit mass is E 2 = \v 2 d . There- 
fore the energy radiated away is E ra d = E\ — E 2 = jf-v 2 d 
per unit mass. The total energy radiated away over the 
age of the bubble is: 

£rad = J ^r 2 cd (p a r P cd v cd )^v 2 d dt (32) 

where the integral is evaluated from t = to t = t. Just 
as we obtained Eq. ^8|, we find that 

-frrad - 7 a x cd g cd 6 Et 



(33) 



Global energy conservation demands that Et = KE^ + 
IEb + KE s h + E rad . With the 4 terms on the RHS scaling 
as a 10+x , the correct normalization of the solution may 
be determined. 

Finally, we can also calculate the fraction of mass lost 
by the star over its lifetime that has yet to pass through 
the inner shock: 



8+ a ( 16 + 3A \ m u 
Psw ~ A/ w ~ a \ 10 + A 



2.4- Solution Procedure 



(34) 



Three parameters fully determine the form of a given sim- 
ilarity solution (A, <j> and x is /x cd ), whilst a further pa- 
rameter, (#), obtained only after a particular solution has 
been found, is a measure of the ratio of mass evaporated 
from the clumps to the mass swept-up from the inter- 
clump medium. For bubbles whose evolution is signifi- 
cantly altered by mass-loading, we require that simulta- 
neously both <j> and 8 are low. 

The similarity equations were integrated with a fifth- 
order accurate adaptive step-size Bulirsch-Stoer method 
using polynomial extrapolation to infinitesimal step size. 
Once the CD was reached, rescaling was implemented with 
the relationships in Eq. ^ so that x cd = 1. The mass, and 
kinetic and thermal energies of the bubble were calculated, 
as were the kinetic energy of the shell and the energy ra- 
diated from it. The correct normalization to satisfy global 
energy conservation was then obtained. Finally, for given 
values of E, Q, and t, x, /, g and h may be used to calcu- 
late the physical variables r, p, u, and e. 
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3. Results 



We first consider similarity solutions obtained by Dyson & 
de Vries ( |1972| ) and Dyson (1973). These earlier works con- 
cerned the dynamical effects of a high velocity stellar wind 
incident on a smooth ambient medium and contain no re- 
sults for mass-loaded flows. Dyson & de Vries originally 
computed the form of the shocked ambient gas behind 
the forward shock by parameterizing the cooling function, 



which demanded that the ambient density, p 



with 



q = — 1. They discovered that the shocked ambient gas 
cooled rapidly, and that this region was very thin com- 
pared with t he overall radius of the bubble. Therefore, 
Dyson ( 1973 ) discarded the shocked ambient region in 
favour of a strong isothermal shock coincident with the 
contact discontinuity, which relieved the restriction on q. 
A central assumption for both papers was that the ve- 
locity of the stellar wind was constant, which led to the 
stellar mass- loss rate varying as M oc t q+2 . Thus for an 
ambient density p oc r~ 2 , the energy input by the stellar 
wind was constant as a function of time. Since in our work 
the stellar wind velocity is constant if the ambient density 
falls off as r~ 2 (i.e. (3 = —2, A = —3), and a constant 
rate of energy input is a central assum ption, we can make 
a direct comparison between Dyson's ( |1973 ) results with 
those obtained here. 

In Fig. H we show results for A = —3 with negligible 
mass- loading (large (f>). We chose specific combinations of 
the variables in order to facilitate comparison with results 
given in Fig. 2 of Dyson ( 1973 ): agreement is excellent. In 
Fig. H we plot the relative radii (or velocities) of the inner 
and outer shocks as a function of g*/g c d-, the ratio of the 
stellar wind velocity to the velocity of the contact disconti- 
nuity (or outer shock) . These calculations were again made 
with negligible mass-loading. The results for the (3 = —2 
case may again be directly compared to the earlier results 
of Dyson (his Fig. 5). In Fig. || we also show results for 
other values of (3. These cannot be compared to Dyson's 
earlier results because when (3^2 the stellar wind veloc- 
ity is no longer constant with time. We find that as the 
value of f3 increases (which also corresponds to the mass- 
loading becoming increasingly dominant at large radii), 
the shocked region broadens for a given ratio of g*/g c d- 

In Fig. |i| we present results for a case in which mass- 
loading dominates the evolution of the bubble. In this par- 
ticular case, A = 4, = 0.2, X[ s = 0.403ir c d, and 9 = 0.62. 
The bubble mass is over 8 times higher than the total 
mass lost by the star, and almost twice the mass of the 
ambient medium swept up into the thin shell. 30 per cent 
of the mass lost by the star has still to pass through the 
inner shock. The high degree of mass-loading increases the 
post-shock density, and dramatically reduces the temper- 
ature in the post-shock region as more and more mass is 
evaporated. This reduces the sound speed, and there is 
eventually a sonic transition close to the CD as the flow 
becomes supersonic in the frame of the clumps. 
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Fig. 3. The ratio of the radii of the inner and outer shocks 
as a function of g*jg c & for cases with negligible mass 
loading. For f3 — —2 the earlier results of Dyson (1973) 
are again recovered (cf. his Fig. 3). For (3^2 the ve- 
locity of the stellar wind is a function of time, so these 
curves cannot be compared to those in Dyson's Fig. 3. 
For (3 = —2,-1,0,1, the corresponding values of A are 
_3 -1 5 4 



Values of several parameters calculated from solutions 
with A = —3 are given in Table @. If the mass-loading 
is negligible (illustrated in the table by 4> = 50,000), we 
find the following: i) as the ratio of Xi S /x c d increases, the 
fraction of energy contained in the form of thermal energy 
in the bubble and that in the form of kinetic energy of the 
shell decreases, whilst that in the form of kinetic energy of 
the bubble increases. The fraction of energy radiated away 
by the shocked interclump gas stays relatively constant, 
and is always a small part (< 6 per cent) of the total 
mechanical energy of the stellar wind over its lifetime, 
ii) The mass contained in the bubble is larger than the 
mass of the swept-up shell only when the relative radius 
of the inner shock is large, iii) For small Xi S /x c d, about 
90 per cent of the total mechanical energy of the stellar 
wind is split approximately evenly between the thermal 
energy contained in the SSW region and the the kinetic 
energy of the swept up shell. As Xi s /x c< i increases, these 
fractions decrease at roughly the same rate, whilst the 
kinetic energy of the gas inside the bubble increases, until 
it eventually may exceed 90 per cent of the total energy, 
iv) The kinetic energy of the swept-up shell is always as 
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Fig. 2. Results for A = —3, Xi s /x c d = 0.359, and 4> = 50, 000. In the limit of large 4> (i- e - negligible mass-loading), the 
earlier results of Dyson (1973) are recovered. Specific combinations of the variables are chosen in order to facilitate 
comparison with results given in Fig. 2 of Dyson (1973). In our work, / is a measure of density, g is a measure of 
velocity, and h is a measure of the thermal energy density. Mach(x) is the Mach number of the flow as a function of x. 



great or greater than the thermal energy of the bubble 
(this is also true for the solutions in Table [ijwith A = — 3, 
0=1.0). 

When we set cj) = 1-0, the bubble is just beginning 
to become mass-loaded for small values of Xi s /x c d- How- 
ever, for A = —3, the swept-up mass dominates the bubble 
mass in this regime, as can be seen from the values in Ta- 
ble [|. We find that solutions for which both <f>b is large 



and Mgh/Afb is small exist only if A > 4. Solutions for 
A = 4 (j3 = 1) and = 0.2 are illustrated in Fig. |. The 
solution with x- m /x c d = 0.420 is strongly mass-loaded, and 
in contrast to the solutions with A = — 3 the mass of the 
swept up shell is also less than Mb. We note that when 
a^is/^cd is varied, notable differences between the profiles 
of /, g, and h occur in Fig. |j| For instance, when the bub- 
ble is only weakly mass-loaded the density at the CD is 
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Fig. 4. Results for a bubble with significant mass-loading. This particular solution is for A = 4, <f> — 0.199, and 
x- ls = 0.403a; c d. The bubble mass is 8.4 times the total mass lost by the star. The ratio of the swept-up mass to the 
bubble mass is M s h/Mb = 0.54. 



smaller than immediately after the inner shock. However, 
once the bubble is more strongly mass-loaded the density 
continuously rises from the inner shock to the CD. Other 
differences include less deceleration of the flow, a sharp 
fall in the temperature (instead of a gradual rise), and a 
rise in the Mach number. Finally, when the bubble be- 
comes strongly mass-loaded, the fractional distribution of 
the total energy between I Eh, KEh, KE s h, and i? ra d for a 
given Xis/xd ratio is substantially different than when the 
mass-loading is negligible. This is hardly surprising given 
the change seen in the profiles of the flow variables. 

An important finding is that for a given A and ratio 
^is/^cdj there exists a minimum value of <p at which we can 



find a solution. In other words, there is a limit to the maxi- 
mum mass-loading which the bubble can undergo. In Fig. |^ 
we show the value of 4>min f° r various values of Xi s /x c d and 
(3. This finding has the corollary that for a given value of 
4> there is a maximum value of a;i s /x c d- The numerical rea- 
son that we are unable to obtain solutions with 4> < 4>min 
is that the denominator of the derivatives from Eqs. 
|l0| approaches zero without the numerators doing likewise. 
The physical reason is that there is a straightforward feed- 
back mechanism limiting the amount of mass-loading that 
may take place: as mass from the clumps is evaporated, 
the temperature of the shocked region falls, which reduces 
the rate at which further mass can be evaporated. Simi- 
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Table 1. Properties of solutions with A = —3 (j3 = —2). We consider two cases: negligible mass-loading (<j> = 50, 000) 
and moderate mass-loading (<p = 1.0). Solutions with A = 4 and <fr = 0.2 are also listed. From left to right, the columns 
indicate: i) the ratio of the shock radii, ii-v) the fraction of energy in the form of thermal energy and kinetic energy in 
the bubble, kinetic energy in the swept-up shell, and radiated energy, vi) the value of 9, vii) the ratio of mass within 
the bubble to the total mass lost by the star, viii) the ratio of mass in the swept-up shell to mass in the bubble, ix) 
the fraction of the mass lost by the star which has yet to pass through the inner shock. 
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larly, if the mass-loading were suddenly inhibited, the tem- 
perature of the shocked region would rise, leading to an 
increase in the evaporation rate. There is no reason to 
believe that actual wind bubbles would be at this limit 
because the degree of mass-loading simply depends on the 
number of clumps present and their individual mass injec- 
tion rates (which in our model is characterized by the value 
of (pb). For a given mass, the number density of clumps, 
n c oc r~ 3 where r c is the clump radius, whilst the mass 
injection rate from each clump rh c oc r c . Thus we can only 
comment that for small clumps, it is more likely that this 
maximum mass-loading limit will actually occur. 

We further note that the radial profile of the Mach 
number of the flow relative to the clumps can take many 
different forms, depending on the input parameters used. 
In Fig. ^ we show several examples. The profile of the 
Mach number is influenced by the fact that mass-loading 
tends to drive the Mach number towards unity, whilst the 
spherical divergence tends to drive it away from unity (cf. 
Hartquist et al. 1986] ). Panel a) is for a case in which mass- 
loading is small. In this case the divergence term wins, and 
drives the Mach number towards zero. Panel b) is for a 
similar situation, but in this case because the inner shock 
is very close to the CD, the initial postshock Mach num- 
ber is greater than unity so the divergence term drives it 
towards infinity. In panel c) we see that the divergence is 
initially dominant, but as the flow approaches the CD the 
Mach number is driven back towards unity. However, this 
behaviour is not due to the flow mass-loading (although 
$b = 2.64), but is rather because the density increases 
near the CD whilst the pressure remains constant, produc- 
ing a decrease in the sound speed (the same effect occurs 



for the flow shown in Fig. |2| which has no mass-loading). 
The same situation occurs in panel d) where the fall in 
sound speed is great enough to drive the flow through 
a sonic transition. The flow corresponding to panel e) is 
very similar to that associated with panel d) except we 
now have two sonic transitions due to the initial post- 
shock flow being supersonic with respect to the clumps. 
In plot f) the flow is significantly mass-loaded (<f>b = 7.86), 
and the mass-loading term dominates the divergence term 
such that the Mach number is driven towards unity. 

Finally, we have calculated profiles of the X-ray emis- 
sivity as a function of radius (see Fig. [§[). We assume that 
the emissivity A oc n 2 ! 1-1 / 2 , which is a good approxima- 
tion over the te mperature range 5 x 10 5 K < T < 5 x 10 7 K 
(cf. Kahn 1976 ). Also plotted in Fig. ||are the radial tem- 
perature profiles. A bubble un- affected by mass- loading 
expanding into a surrounding medium with an r~ 2 profile 
of density has a higher central temperature but a lower 
central emissivity than at its limb. Conversely, if it is ex- 
panding into a medium with an r 1 profile of density, the 
situation is reversed, as the central temperature is lower 
and the central emissivity is higher than at its limb. If the 
bubble is mass-loaded, the general trend is for a reduction 
in the central emissivity and an increase in the central 
temperature relative to the limb. To date, the only stel- 
lar wind bubbles which have been successfully observed in 
X-rays are two Wolf-Rayet ring nebulae, NGC 6888 (e.g. 
Wrigge et al.|1998[) and S308 (Wrigge |1999[). Both have 



X-ray luminosities lower than expected from the standard 



model (Weaver et al. 1977). This discrepancy is normally 



attributed to two possibilities : either conductive evapora- 
tion of gas from the cold dense outer shell into the bubble 
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Fig. 5. Results for bubbles in which the radial dependence of the mass- loading is oc r 4 . Three solutions with different 
values of Xi S /x c d are shown: 0.177, 0.374, and 0.420 (solid, dotted, and dashed lines). <&b is equal to 1.26, 1.36 and 
7.86 respectively, whilst M^/M^ is equal to 4151, 171, 0.71 respectively. For ease of comparison, / is scaled such that 
its pre-shock value is unity in each case. 



interior, or enhanced cooling resulting from high metal- 
licities in the cooling gas (see MacLow 2000 ). It is not 
surprising, therefore, that our mass-loading simulations 
indicate that evaporation from clumps within the bubble 
may also be compatible with current observations. The 
latest X-ray satellites, Chandra and XMM, should provide 
significantly improved observations which may enable us 
to distinguish between these competing mechanisms. One 
might also expect that they may allow discrimination be- 
tween a mass-loaded bubble expanding into an ambient 
medium with par 1 and a bubble unaffected by mass- 
loading expanding into a medium with p cx r~ 2 . 



4. Summary 

We have investigated the evolution of a mass-loaded wind 
blown bubble with a constant rate of ejection of mechani- 
cal energy from a central stellar wind. Mass-loading is as- 
sumed to occur by the conductively induced evaporation 
of clumps in the region of shocked stellar wind. A cen- 
tral assumption is that the shocked interclump medium 
rapidly cools to a dense, negligibly thin, shell through 
which the clumps can pass undisturbed. The requirement 
that the solution be self-similar imposes a link between 
the radial variation of the interclump density (p oc r@) 
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Fig. 6. The minimum value of <p (i-e- corresponding to the maximum possible mass-loading) for which solutions could 
be found, as a function of Xi S /x c d and (3 (or, alternatively A). 



and the mass- loading from the clumps (p oc r A ) which 
forces A= (5 + 7/3)/3. 

We first reproduced solutions with negligible mass- 
loading and (3 = —2, which we compared with results 



obtained by Dyson (1973). Excellent agreement for the 
structure of the bubble was found. We also found that the 
ratio of the radii (velocities) of the inner and outer shocks 
decreases with increasing (3 for a given ratio of the current 
stellar wind velocity to the CD velocity. 

We then investigated the changes in the structure of 
the bubble for different values of A, (f>, and Xi s /x c( \- The 
central conclusions are: 



— Substantial mass-loading of the wind-blown bubble 
can occur over a wide range of A. However, to addi- 
tionally satisfy the requirement that the bubble mass 
is larger than the mass of the swept-up shell, large val- 
ues of A are needed (A > 4). 

— The profiles of the flow variables are significantly al- 
tered under conditions of large mass-loading. In par- 
ticular, the average density of the shocked region is 
larger, the deceleration of the flow is shallower, and the 
temperature of the shocked wind rapidly decreases. 

— The Mach number of the shocked stellar wind relative 
to the clumps which are injecting the mass can take 
several different forms. The shocked region can be ei- 



Pittard, Dyson & Hartquist: Self-similar evolution of Mass-Loaded WBBs 



11 





0.940 0.950 0.960 0.970 0.980 0.990 1.000 





'I I I I 1 



(0 



0.90 [ i i i i I i i i i I i i i i I i i i i I i i i i 1 0.50 n I 

0.70 0.75 0.80 0.85 0.90 0.95 1.00 0.75 0.80 0.85 0.90 0.95 1.00 0.40 0.50 



"I I I 



0.60 0.70 0.80 0.90 



1.00 



Fig. 7. The Mach number as a function of x for a number of different solutions. With respect to the clumps, the 
shocked region can be either entirely subsonic, entirely supersonic, or have up to two sonic points. The corresponding 
solution shown in each panel is: a A = 4, <j> = 0.2, Xi s /x c< i = 0.177; b A = —3, <fi = 1.0, Xi s /x c d = 0.950; c A = —3, 
<j> = 1.0, x is /x cd = 0.120; d A = -3, 4> = 50,000, a; is /x cd = 0.688; e A = -3, 4> = 50,000, x is /x cd = 1.301; f A = 4, 
4> = 0.2, x is /x cd = 0.420. 



ther entirely subsonic, entirely supersonic, or have one 
or sometimes two sonic points. 

— For a given A and ratio of the shock radii, there exists 
a maximum mass- loading that can occur. 

— Mass-loading tends to reduce the emissivity in the in- 
terior of the bubble relative to its limb, whilst simulta- 
neously increasing the central temperature relative to 
the limb temperature. 
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